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Introduction 


Circular waveguide-fed apertures or conical horns are often 
used as elements of planar array antennas. The purpose of this 
present paper is to document a computer code which calculates the 
mutual coupling between circular apertures in a conductive plane. 
The program is quite general in that the apertures do not have to 
be the same sizes, nor do they have to be polarized in the same 
direction. In addition, several waveguide modes (TE and/or TM) 
may be specified in the apertures and the mutual coupling between 
all combinations of apertures and modes will be calculated. The 
program also allows multiple layers of homogeneous dielectrics to 
be placed over the aperture array. Outside the layered region, 
one can specify either a homogeneous half-space or a perfect 
reflecting surface. The program is a revision of an earlier one 
(ref.l) in an attempt to overcome some difficulties encountered 
with the earlier version and to expand its versatility. 
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Symbols 


E J 

J (z) 


J'(z) 

m 


N+l 


m 


m 


n 


n 


J 

N 

R 

R||(0) 

Rj_(0) 

RnO 3 ) 


RiO) 

s 

ij 

[S] 

TE 

TM 

x,y,z 


Radius of i aperture. 

Radius of j aperture. 

th 

Thickness of n homogeneous dielectric layer. 

i aperture electric field polarization vector. 

j aperture electric field polarization vector. 

Bessel function of the first kind of order m and 
argument z . 

Derivative of J (z) with respect to z. 

m 

v=i 

Wave propagation constant in free space, 2n/X. 

Wave propagation constant in n layer. 

Wave propagation constant in homogeneous half-space. 
First index of waveguide mode in i aperture. 

th 

First index of waveguide mode in j aperture. 

Second index of waveguide mode in i aperture. 

Second index of waveguide mode in j aperture. 

Number of homogeneous layers outside array plane. 

Center-to-center spacing between apertures . 

Reflection coefficient for plane wave with electric 
field polarized parallel to plane of incidence. 

Reflection coefficient for plane wave with electric 
field polarized perpendicular to plane of incidence. 

Plane wave reflection coefficient for parallel 
polarization for exterior scattering problem. 

Plane wave reflection coefficient for perpendicular 
polarization for exterior scattering problem. 

Mutual coupling between i and j apertures for 
one mode in each aperture th 

Complex square matrix whose elements are S tj . 

Denotes transverse (to z) electric field. 

Denotes transverse (to z) magnetic field. 

Cartesian coordinate variables. 
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x-coordinate of i aperture center. 

th 

x-coordinate of 1 aperture center. 

th 

y-coordinate of i aperture center. 

y-coordinate of j aperture center. 

Mutual admittance between i and j fch apertures 
for one mode in each aperture. 1 

Complex square matrix whose elements are Y tj . 

Complex diagonal matrix whose elements are the 
waveguide modal characteristic admittances. 

Normalized radial propagation constant in 
cylindrical coordinate system. 

Distance from outer surface of N t dielectric layer 
to exterior scattering boundary. 

Permittivity of free space. 

Permittivity of region just outside aperture plane. 

Permittivity of n homogeneous layer. 

Permittivity of homogeneous half-space outside of 
N th la y er * 

Equivalent permittivity just inside of exterior 
scattering boundary. 

Wavelength in free space. 

Permeability of free space. 

Permeability of region just outside of aperture. 

Permeability of n homogeneous layer. 

Permeability of homogeneous half-space outside of 
N layer. 

t h 

Equivalent permeability just inside of exterior 
scattering boundary. 

Polarization of i aperture with respect to y-axis. 

t h 

Polarization of j aperture with respect to y-axis . 

t h 

V V 

n zero of J (z). 

th m v ’ 

n zero of J'(z). 

th m v ' 
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Theory 


The analytical formulation for the electromagnetic 
interaction between apertures in a planar array has been developed 
earlier (ref.l) and will not be repeated here. 

The mutual interaction between the various modes among all 
the apertures in the array is obtained from the appropriate 
coefficients of the scattering matrix, [S] , as follows: 

[S] - [ [Y.] - [Y] ][ [Y.] ♦ [Y] j- 1 (1) 

Where [Yo] is a diagonal matrix whose elements are the 
characteristic admittances of the waveguide modes and [y] is the 
admittance matrix whose elements are the mutual admittances 
between each waveguide mode and all modes in each waveguide 
aperture. Therefore, for example, if the array consisted of 5 
apertures and each aperture field was assumed to be the 
superposition of 7 waveguide modes (TE and/or TM) , then (1) would 
consist of square complex matrices each of size 35 by 35. The 
sizes of the matrices can expand very rapidly if one does not use 
some "engineering judgment" in selecting the modes to approximate 
the aperture fields. The expressions for calculating the elements 
of the admittance matrix (i.e., the complex mutual admittance 
between pairs of waveguide modes ) are derived in reference 1 . 

Y „ - -/W /V-j I" { *,<*> »„«*> 

- W 2 (S) C,(0) <,((*) V U (P) | e dp (2) 

where , 

7=1 for m = 0 

m i 1 

= 2 otherwise. 
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W i (/3) and W 2 (/9) are defined as 


w a O) 


1 - R.(P) 

1 + R„(j3) 



1- 



w 2 0) 


1 - R ± (f?) 

1 + RjW 



(3) 

(4) 


where (for real angles) /3 can be interpreted as the sine of the 
angle of incidence for a plane wave with R| ( (/3) and Rj_(/3) being the 
reflection coefficients in the outward direction evaluated at the 
aperture boundary for parallel and perpendicular polarization at 
an air-dielectric interface. Notice from equation 2 that the 
mutual admittance is a weighted sum over all angles of incidence, 
including real space (O^/Ssl) and invisible space (j3>l). An 
alternate form for W j (/3) and W z (£) is 


W^/3) 


w 2 «3) 


f ’ jk o C i ) f 

l C o / [ 

( i f f ;^>] 

t“=3Wi 


(5) 


( 6 ) 


where g^P) and f j (j3) are the solutions to the wave equations in 
the external region evaluated at the aperture plane, and g'(/3) and 
f'(|9) are the corresponding derivatives in the direction normal to 
the aperture plane. c i and are the permittivity and 
permeability of the external medium evaluated just outside the 
aperture plane. 

By applying the boundary conditions of continuity of 
tangential fields at the interface of each homogeneous dielectric 
layer, one can evaluate the mutual admittance between aperture 
modal fields with a multilayer dielectric medium outside (as shown 
in figure 1 ) . 
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Figure la. Coordinate geometry for the i^ h and j fch elements of 
a planar array of circular waveguide-fed apertures 


Figure lb. Cross-section of i and j circular waveguide-fed 
° th th 

apertures radiating into a multi-layer dielectric 
with an exterior arbitrary reflecting boundary. 




* 


Applying boundary conditions at layer interfaces, results in the 
following recursive equations. 



for n = N, (N-l ), (N-2 ),..., 3,2 , 1 . where is defined as 



i e * k/k o 
; e > k/k o 


Then starting with 




9; +1 o)l 

v 



! V S k , 
• V » 


( 9 ) 


equations 7 and 8 are evaluated successively from the outermost 
homogeneous layer to the layer just outside of the aperture plane. 
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In the special case of a perfect electric conductor of 
infinite extent on the outer surface of the N layer, f (/9)=0 
and g' O)=0; therefore, the starting conditions for equations 7 

N+l 

and 8 become 


ff^r 

— — ]r 4 

cos(k d )' 

z N 

N j 

” * A 4 

Z 

sin <’VVJ 


r 1C ^ 

f Sln(k z d » ) ] 


A 

Z 

[ cos(k 2 d N ) 


(10a) 


(10b) 


If the medium outside of the N dielectric layer is not a 

t h 

homogeneous half-space, the starting conditions (equation 9) are 
replaced by 



= -j K((3) 

M I 

P 

r 1 - R£(/3) exp(-j26 K(|3)) 

J 

ll 

N+l 

V 

k 1 + R±(/3) exp(-j26 K(j3)y 


- -3 K(/3) 

c 

p 

' 1 + Rj(/3) exp(-j26 K(p))' 

\ > 

C N+1 

v. 

k 1 - RnO) exp(-j25 K(|3)) > 


(Ha) 


(lib) 


where. 


K(/3) 


/ ^ 

-J 


; k /3 s k 

' QT H+l 


; k £ > k 
' cr h+i 


r£(/ 3) and R^O) are the plane wave reflection coefficients, 
evaluated at the external scattering boundary (shown dashed in 
figure lb) for perpendicular and parallel polarization of the 
electric field vector with respect to the plane of incidence, c 

p 

and M p are the equivalent permittivity and equivalent permeability 
just inside of the external scattering boundary. The starting 
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conditions in equation 11 allow for future modifications of the 
computer code to include problems for which the plane wave 
reflection is known from solutions to exterior scattering 
problems, such as objects in the vicinity of the array or an 
inhomogeneous ionized plasma. 

The other quantities in equation 2 are related to the Fourier 
transforms of the aperture modal fields as derived in reference 1. 

For the TE modes in the i aperture: 


(*; . ) ! (v.) j ; (v/) 

V 1 i J l 

{(*-,“,) - ( k 0 a / } ^ (*%",) 



03) 


( k o a ,) m i MW)/ Ow) 



(13) 


For the TM inodes in the i aperture: 

<™«J) * 0 


e™(f) = 


( k o a .) ( k oV) J . (V.f) 

, 1 

~ ( k o a / 


(14) 

(15) 


For the j aperture, the i subscripts are replaced with j in 
equations 12-15. 

If R=0, then the i^and j apertures are coincident and 
equation 2 is the mutual admittance between modal fields in the 
same aperture. When R=0, the orthogonality of the modal functions 
results in U 1 j (/3)=V ij O)= 0 except when m^ny When R=0 and m^m ; 
U tJ 03) and V tj (/3) are as given in Table I. 
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Table II 
( R > 0 ) 


Modes 

u,,(P) 

v, ,(e) 

TE t TEj 

(-1) ' ( C + B + (0) 
m. 

- (-1) c_ B_(f3) ) 

m i 

(-1) J ( C + B + (P) 
m. 

+ (-1) C_ B J/3) ) 

TM t TMj 

-(-1) J ( C + B + (0) 
m 

+ (-1) C_ B_«3) ) 

0 

TE j TMj 

( S + B + (0) 
m 

" (-1) S. B_(/3) ) 

0 

TM t TEj 

(-1) 1 ( S + B + (B) 

m ! 

- (“1) S_ B_(S) ) 

0 
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Appendix 

Computer Program 


The computer code is written in FORTRAN IV and uses double 
precision variables throughout, except for integers. It was found 
that double precision was necessary in order to obtain sufficient 
accuracy on a VAX computer. The input data is read from logical 
unit 5. Results are output to logical unit 6. A flow chart 
indicating the subroutine calls is given and also a description of 
the input parameters with sample input and output data. The 
computer code has been verified against results in reference 1. 


Subroutine Flowchart 



A-l 





















Definition of Input Data 


NHOLE 

NMODE 

NUMTE 

NUMTM 


MIJ(I) 

NIJ(I) 

MIJP(I) 

NIJP(I) 

AIJ(I) 

XIJ(I) 

YIJ(I) 

PHIJP(I) 

F 

FCON 

ER 

CEP 

CUP 

NLAY 


D( I ) 
CE( I ) 
CU(I) 


Number of apertures in array. 

Number of modes per aperture. 

Number of TE-modes per aperture. 

Number of TM-modes per aperture. 

(NOTE: same modes assumed for all apertures). 

First index of I-th TE mode. 

Second index of I-th TE mode. 

First index of I-th TM mode. 

Second index of I-th TM mode. 

Radius of I-th aperture, 
x-coordinate of center of I-th aperture, 
y-coordinate of center of I-th aperture. 
Polarization of I-th aperture (degrees CCW). 

Frequency (Hertz). 

Factor for converting input to centimeters. 
Dielectric constant of waveguide interior. 

Relative epsilon outside of layered region. 
Relative mu outside of layered region. 

Number of homogeneous dielectric layers. 

Thickness of I-th homogeneous layer. 

Relative epsilon of I-th layer. 

Relative mu of I-th layer. 


Input Data 


(1) Read NHOLE, NMODE, NUMTE, NUMTM 
Check for end-of-file on unit 5. 

(2) If NUMTE=0 , skip to (3). 

Read MIJ(I),NIJ(I) ; 1=1 To NUMTE 

(3) If NUMTM=0, skip to (4). 

Read MI JP( I ) , NI JP( I ) ; 1=1 To NUMTM 

(4) Read AIJ ( I ) , XI J ( I ) , YI J ( I ) , PHI JP( I ) ; 1=1 To NHOLE 

(5) Read F, FCON, ER 

(6) Read CEP, CUP 

(7) Read NLAY 

(8) If NLAY=0, skip to (9). 

Read D(I),CE(I),CU(I) ; 1=1 To NLAY 

(9) Perform calculations. 

Read next set of input data (1). 
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Sample Input Data 


2 110 
11 

0.75 0.0 0.0 0.0 
0.75 0.0 2.5 0.0 
6.0D09 2.54 1.0 
( 1 . 0 , 0 . 0 ) ( 1 . 0 , 0 . 0 ) 

1 

0.18 (2.6,-0.0156) (1.0, 0.0) 


Sample Output Data 


MUTUAL COUPLING OF CIRCULAR APERTURES 
RADIATING INTO MULTI-LAYERS 
NUMBER OF LAYERS = 1 


LAYER THICKNESS 
1 0. 1800 


COMPLEX EPSILON 
2.60000, -0.01560) 


COMPLEX MU 
1 . 00 , 0 . 00 ) 


OUTSIDE HALF-SPACE EPSILON - 1.00000 0.00000 

OUTSIDE HALF-SPACE MU = 1.00000 0.00000 


**** INPUT DIMENSIONS IN INCHES **** 
*»** FREQUENCY = .600000E+10 HERTZ **** 


++++++++ APERTURE ARRAY GEOMETRY +++++++++ 


HOLE RADIUS X 

1 0. 7500 0. 000 

2 0. 7500 0. 000 


Y POL. 

0. 000 0. 000 

2.500 0.000 


MODE 1 = TE-11 


CY( 1, 1)=( 0. 3415E-02, 0. 1691E-02) BB= 17.00 
CY( 1, 2)=( 0. 3443E-04, -0. 3158E-03) BB= 15.00 


IH= 1 IM= 1 JH= 1 JM= 1 
IH= 1 IM= 1 JH= 2 JM= 1 


YMN( 1)=( 0. 1695E-02, 0.0000E+00) IHOLE= 1 IMODE= 1 

YMN( 2)=( 0. 1695E-02, 0.0000E+00) IHOLE= 2 IM0DE= 1 


SCATTERING MATRIX 
S( 1, l)*(-0.4036E+00,-0. 1964E+00) 
S( 1, 2)=( 0. 1871E-01, 0. 3199E-01 ) 
S( 2, 1)=( 0. 1871E-01, 0.3199E-01) 
S( 2, 2 ) = ( — 0. 4036E+00, -0. 1964E+00) 


-6.9570 DB 
-28.6224 DB 
-28.6224 DB 
-6.9570 DB 


-154.0525 DEG 
59.6795 DEG 
59.6795 DEG 
-154. 0525 DEG 
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£*******•********* £yQ It********************************************** 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


Mutual admittance and scattering matrix calculation for 
circular waveguide array with TE-MN and TM-MN modes. 

NOTE: Field equations of Marcuvitz show 

TM modes polarized 90 degrees to TE modes, 
program modified to rotate polarization of TM modes 


by PI/2 such that TE-11 and TM-11 are polarized the same.* 


Program can handle multiple homogeneous layers over apertures. 
(INPUT on unit #5) (OUTPUT on unit #6) 


PROGRAM: 


M.C. Bailey 
Hampton, VA 


(1989) * 


(-I********************************************************************** 


102 FORMAT (IX, 79 (1H+)) 

WRITE(6, 102) 

CALL CONSTS 
100 WRITE(6, 102) 

CALL INPUT 
CALL YMAT 
CALL SMAT 
GO TO 100 
END 
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SUBROUTINE CONSTS 

IMPLICIT REAL*8 (A, B, D-H, O-Z) 

IMPLICIT COMPLEX* 16 (C) 

COMMON /ZEROS/ XMNP (9, 6) , XMN (9, 6) 

COMMON /CONST/ CJ, PI , TWOPI , DTOR, RTOD, BIG, SMALL, BMAX, YO, MAX, LAYM 
DATA XMNP(1, l),XMNP(l,2),XMNP(l,3)/3. 831706,7. 0155867, 10. 173468/ 
DATA XMNP(2, 1 ) , XMNP(2, 2) , XMNP(2, 3)/l. 84118, 5. 33144, 8. 53632/ 

DATA XMNP(3, 1 ) , XMNP(3, 2) , XMNP(3, 3)/3. 05424, 6. 70613, 9. 96947/ 

DATA XMNP(4, 1 ) , XMNP(4, 2) , XMNP(4, 3)/4. 20119, 8. 01524, 11. 34592/ 

DATA XMNP(5, 1 ) , XMNP(5, 2) , XMNP(5, 3)/5. 31755, 9. 28240, 12.68191/ 

DATA XMNP(6, 1 ) , XMNP(6, 2) , XMNP(6, 3)/6. 41562, 10.51986, 13.98719/ 
DATA XMNP(7, 1 ) , XMNP(7, 2) , XMNP(7, 3)/7. 50127, 11.73494, 15.26818/ 
DATA XMNP(8, 1 ) , XMNP(8, 2) , XMNP(8, 3)/8. 57784, 12.93239, 16.52937/ 
DATA XMNP(9, 1 ) , XMNP(9, 2) , XMNP(9, 3)/9. 64742, 14. 11552, 17.77401/ 
DATA XMNP( 1 , 4) , XMNP( 1 , 5) , XMNP( 1 , 6 )/13. 323692, 16.47063, 19.615859/ 
DATA XMNP(2, 4) , XMNP(2, 5) , XMNP(2, 6)/ll. 70600, 14.86359, 18.01553/ 
DATA XMNP(3, 4) , XMNP(3, 5) , XMNP(3, 6)/13. 17037, 16.34752, 19.51291/ 
DATA XMNP (4,4), XMNP (4,5), XMNP (4,6)/14.58585,17. 78875 , 20 . 97248/ 
DATA XMNP (5, 4) , XMNP (5, 5) , XMNP (5, 6)/15. 96711, 19. 19603,22. 40103/ 
DATA XMNP (6,4), XMNP (6, 5) , XMNP(6, 6)/17. 31284, 20. 57551, 23. 80358/ 
DATA XMNP (7, 4) , XMNP (7, 5) , XMNP (7, 6)/18. 63744, 21 . 93172, 25. 18393/ 
DATA XMNP ( 8 , 4) , XMNP(8, 5) , XMNP(8, 6)/19. 94185, 23. 26805, 26. 54503/ 
DATA XMNP (9,4), XMNP (9,5), XMNP (9, 6)/21. 22906, 24. 58720, 27. 88927/ 
DATA XMN( 1, 1 ) , XMN( 1 , 2) , XMN( 1, 3)/2. 4048256,5. 5200781,8.6537279/ 
DATA XMN(2, 1),XMN(2, 2), XMN(2,3)/3. 8317060, 7. 0155867, 10. 1734681/ 
DATA XMN(3, 1 ) , XMN(3, 2) , XMN(3, 3)/5. 1356223,8.4172441, 11.6198412/ 
DATA XMN(4, 1 ) , XMN(4, 2) , XMN(4, 3)/6. 3801619, 9. 7610231, 13. 0152007/ 
DATA XMN ( 5 , 1 ) , XMN ( 5 , 2 ) , XMN ( 5 , 3 ) /7 . 5883427 , 1 1 . 0647095 , 1 4 . 3725367/ 
DATA XMN (6, 1 ) , XMN(6, 2) , XMN(6, 3)/8. 7714838, 12. 3386042, 15. 7001741/ 
DATA XMN(7, 1),XMN(7, 2), XMN(7,3)/9. 93611, 13.58929, 17.00382/ 

DATA XMN (8, 1 ) , XMN (8, 2) , XMN (8, 3)/l 1 . 08637, 14.82127, 18.28758/ 

DATA XMN (9, 1 ) , XMN (9, 2) , XMN (9, 3)/12. 22509, 16.03777, 19.55454/ 

DATA XMN (1,4), XMN (1,5), XMN ( 1 , 6)/l 1 . 7915344, 14.9309177, 18.071064/ 
DATA XMN (2, 4) , XMN (2, 5) , XMN (2, 6)/13. 32369, 16.47063, 19.61586/ 

DATA XMN (3, 4) , XMN (3, 5) , XMN (3, 6 )/14. 79595, 17.95982,21. 11700/ 

DATA XMN ( 4 , 4 ) , XMN ( 4 , 5 ) , XMN ( 4 , 6 ) /I 6 . 22347 , 1 9 . 40942 , 22 . 58273/ 

DATA XMN(5,4),XMN(5,5),XMN(5,6)/17. 61597, 20. 82693, 24. 01902/ 

DATA XMN(6, 4) , XMN(6, 5) , XMN(6, 6)/18. 98013, 22. 21780, 25. 43034/ 

DATA XMN ( 7 , 4 ) , XMN ( 7 , 5 ) , XMN ( 7 , 6 ) /20 . 32079 , 23 . 58608 , 26 . 820 1 5/ 

DATA XMN (8,4), XMN (8, 5) , XMN(8, 6 )/21 . 64154, 24. 93493, 28. 19119/ 

DATA XMN (9, 4) , XMN (9, 5) , XMN (9, 6) /22. 94517, 26. 26681, 29. 54566/ 

DATA CJ, BIG, SMALL/ (O.ODO, 1 . ODO) , 1 . OD+38, 1.0D-38/ 

DATA BMAX, MAX, LAYM/200. 0, 40, 600/ 

PI=2. 0*DASIN ( 1 . ODO ) 

Y0=1. 0/(120. 0*PI) 

DTOR=P 1/180.0 
RTOD= 180. O/PI 
TWOPI =2. 0*PI 
RETURN 
END 
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SUBROUTINE INPUT 


c: Reads input data 


IMPLICIT REAL* 8 (A, B, D-H, O-Z) 

IMPLICIT COMPLEX* 16 (C) 

LOGICAL LPOL, LSIZE, LSAME, LCON, LPLAZ 

COMMON DIA, MMM, NHOLE, NMODE, NUMTE, NUMTM, F, FCON, ER, LPOL, LSIZE, LSAME 
COMMON /LAYER/ CEP, CUP, NLAY, D(601 ), CE(601 ), CU(601 ), CEUP, LCON, LPLAZ 
COMMON /CONST/ CJ, PI , TWOPI, DTOR, RTOD, BIG, SMALL, BMAX, YO, MAX, LAYM 
COMMON /ARRYS/ MIJ(20) , NI J(20) , MIJP(20) , NI JP(20) , 

A CY(40, 40) , CA(40, 40) , CB(40, 40) , IPIV(40) , INDX(40, 2) , 

B AIJ(40) , XIJ(40), YIJ(40), PHIJP(40) 

101 FORMAT! IX) 

102 FORMAT ( IX, 79 (1H+)) 

Q****f*t**f*t*»**»**t***» INPUT DATA ********************************* 

C************************************************************* ********** 


C NHOLE = Number of apertures in array. * 
C NMODE = Number of modes per aperture. * 
C NUMTE * Number of TE-modes per aperture. * 
C NUMTM = Number of TM-modes per aperture. * 
C NOTE: ( NUMTE+NUMTM ) =NMODE * 
C The same modes are assumed in all apertures. * 


READ ( 5 , * , END=9999 ) NHOLE , NMODE , NUMTE , NUMTM 

MMM=NHOLE*NMODE 

IF(MMM.LE.MAX) GO TO 120 

WRITE(6, 110) 

110 FORMAT! IX’ (NHOLE* NMODE) EXCEEDS DIMENSION OF CY’ ) 

STOP 

120 WRITE (6, 130) 

130 FORMAT! IX’ MUTUAL COUPLING OF CIRCULAR APERTURES’ ) 

IF(NUMTE.GT. 20) GO TO 500 
IF (NUMTM. GT. 20) GO TO 500 
IF(NUMTE.GT. NMODE) GO TO 500 
IF(NUMTM.GT. NMODE) GO TO 500 
IF(NM0DE. GT. 40) GO TO 500 
IF! (NUMTE+NUMTM). NE. NMODE) GO TO 500 
IF (NUMTE. EQ. 0) GO TO 150 

£*************»*********** ********************************** ************ 


C MIJ ( I ) , NI J ( I ) = M,N indices of I-th TE mode. * 

READ (5, 140) ( (MI J( I ) , NI J( I ) ) , 1-1, NUMTE) 

140 FORMAT (20(211, IX)) 

150 IF(NUMTM.EQ.O) GO TO 170 

£*•*•************•**********•*•**************•*************«******»**,',** 


C MIJP( I ) , NIJP( I ) = M, N indices of I-th TM mode. * 

READ (5, 160) ( (MI JP ( I ) , NIJP( I ) ) , 1=1, NUMTM) 

160 FORMAT(20 (21 1 , IX) ) 
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170 CONTINUE 


******************* 


C 

c 

c 

c 


AIJ(I) 

XIJ(I) 

YIJ(I) 

PHIJP(I) 


Radius of I-th aperture. * 
X-coordinate of center of I-th aperutre. * 
Y-coordinate of center of I-th aperture. * 
Angular rotation of Xl-axis WRT X-axis (degrees CCW). * 


READ (5, * ) ( ( AI J ( I ) , XI J ( I ) , YI J ( I ) , PHI JP( I ) ) , 1=1 , NHOLE ) 

«*»*** **»»*****»»*»*»**«****»***«•**»***»*»•***»»***»**** *»**»***»» 

C F = Frequency (Hertz). * 

C FCON = Factor for converting input dimensions to centimeters. * 

C (2.54 for input in inches) (1.0 for input in centimeters).* 

C ER = Relative dielectric constant of material filling waveguides 


C 

C 

c 


For input dimensions in wavelengths set: 
F=3. 0E10 
FC0N=1. 0 


* 

* 

* 

u* 


READ(5, *)F, FCON, ER 

£******»****«********»********»*»**************************»****»******* 


c 

c 

c 

c 

c 

c 

c 


CEP = Complex (relative) Epsilon of half-space outside of layers.* 
CUP = Complex (relative) Mu of half-space outside of layers. * 


Input format is: (A, -B) (C, -D) 

For free-space outside: (1.0, 0.0) (1.0, 0.0) 

For a metal sheet approximation: (-1. 0E10, 0. 0) (1.0, 0.0) 

or CABS ( CEP*CUP ) greater than 1.0E09 


READ (5, * )CEP, CUP 

£***»****•**•***********«**«****•**»************»*•*•*******•«******«**« 


NLAY = Number of homogeneous layers over apertures. 


NLAY=0 


(half-space) . 


READ (5,*) NLAY 

IF(NLAY. GT. LAYM)ST0P 4444 

IF(NLAY.EQ.O) GO TO 210 


Ml******************** 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


NOTE: If NLAY=0, do not inter the following data. 

D ( I ) = Thickness of I-th layer. 

CE(I) = Complex Epsilon of I-th layer (relative to free space). 
CU(I) = Complex Mu of I-th layer (relative to free space). 

NOTE: Values of CE(I) must be input with format: (A, -B) 

Values of CU(I) must be input with format: (C,-D) 

where dielectric loss tangent = B/A 


* 

* 

* 

* 

* 

* 

★ 

* 

* 

* 
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C CAUTION: * 

C An input value of 0.0 for B could produce erroneous results. * 

C A value of 0. 0001 for the dielectric loss tangent should yield * 

C results which approximate a lossless dielectric in most cases. * 

C (Suggest varying B to determine sensitivity to small values) * 

DO 180 1=1, NLAY 

READ (5,*)D(I),CE(I),CU(I) 

180 CONTINUE 

WRITE (6, 190) NLAY 

190 FORMAT (IX’ RADIATING INTO MULTI-LAYERS’ /IX’ NUMBER OF LAYERS = ’14/) 
WRITE(6, 200) 

200 FORMAT (IX’ LAYER ’ 3X ’ TH I CKNESS ’ 5X ’ COMPLEX EPS I LON ’ 1 OX ’ COMPLEX MU ’ ) 
210 IF(DABS(ER) . LT. 1. 0D-04)ER=1. 0 

IF(DABS(FC0N) . LT. 1. OD-04)FCON=1. 0 
XLAM= ( 30 . /FC0N ) / ( F* 1 . 0D-09 ) 

IF(DABS(XLAM-1. 0) . LT. 0. 0001 )XLAM=1. 0 

OMEGA=TWOP I *F 

IF (NLAY. EQ. 0) GO TO 280 

DO 270 1=1, NLAY 

WRITE (6, 220)1 ,D(I),CE(I),CU(I) 

220 F0RMAT( 1X14, 2XF9. 4, 3X’ (’F10. 5’ , ’F9. 5’ )’3X’ (’F6.2’ , ’F6. 2’ )’ ) 

D ( I ) =D ( I ) /XLAM 
EPR=DREAL ( CE ( I ) ) 

EPI=DIMAG(CE( I ) ) 

IF(EPI. LE. -0. 0001*EPR)GO TO 270 
IF(EPI. LT. 0. 0)GO TO 240 
EPI=-EPI 
WRITE(6, 230)1 

230 FORMAT (IX’ NOTE: The imaginary part of the dielectric constant’/ 

A8X’for the layered region must be negative.’/ 

B8X’This has been corrected for layer ’12/) 

240 IF(EPI. LE. -0. 0001*EPR)GO TO 260 
TANL=DABS ( EPI/EPR ) 

WRITE (6, 250) I, TANL 

250 FORMATf IX’ CAUTION: The dielectric loss tangent for the layered’/ 

A11X’ region must not be zero. Erroneous results could also be’/ 
B11X’ obtained for very small values of dielectric loss.’/ 

CllX’The loss tangent for layer ’12’ is 'Ell. 4/) 

260 CE( I )=DCMPLX(EPR, EPI ) 

270 CONTINUE 

WRITE(6, 101) 

280 CONTINUE 
290 LCON=. FALSE. 

CEUP=CEP*CUP 

IF(CDABS(CEUP).GE. 1.0D+9) LC0N=. TRUE. 

ERC=DREAL ( CEP ) 

EIC=DIMAG(CEP) 

NSIGN=1 

IF(ERC. LT. 0. 0)NSIGN=-1 
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ERCP=NSIGN*ERC 

IF(ERCP. GT. 1. 0D+10)ERCP=1 . 0D+10 

ERC=NSIGN*ERCP 

CEP=DCMPLX (ERC, EIC) 

IFCNLAY.EQ.O. AND. ( . NOT. LCON) . AND. (. NOT. LPLAZ) )WRITE(6, 300) 
300 FORMAT ( IX’ RADIATING INTO HALF-SPACE’/) 

IF(LC0N)WRITE(6, 320) 

310 CONTINUE 

320 FORMAT! IX’**** PERFECT CONDUCTOR OUTSIDE ****’/) 

I F ( . NOT . LCON ) WR I TE ( 6 , 330 ) CEP , CUP 
330 FORMAT! IX’ OUTSIDE HALF-SPACE EPSILON = ’2F18.5/ 

+ IX’ OUTSIDE HALF-SPACE MU = ’2F18.5/) 

IF(FCON. EQ. 2. 54)WRITE(6, 340) 

340 FORMAT! IX’**** INPUT DIMENSIONS IN INCHES ****’) 

IF ( (FCON. EQ. 1.0). AND. (F. EQ. 3. 0E10) ) WRITE (6, 350) 

350 FORMAT! IX’**** INPUT DIMENSIONS IN WAVELENGTHS ****’) 

IF! (FCON. EQ. 1.0). AND. (F. NE. 3. 0E10) )WRITE(6, 360) 

360 FORMAT! IX’**** INPUT DIMENSIONS IN CENTIMETERS ****’) 

IF! (FCON. NE. 1. 0) . AND. (FCON. NE. 2. 54) )WRITE(6, 370)FCON 
370 FORMAT! IX’**** INPUT CONVERSION FACTOR = ’F10.5’ ****’) 

IF! (FCON. NE. 1.0). OR. (F. NE. 3. 0E10) )WRITE(6, 380)F 
380 FORMAT! IX’ **** FREQUENCY =’ Ell. 6’ HERTZ ****’) 

WRITE(6, 390) 

390 FORMAT! /IX’ ++++++++ APERTURE ARRAY GEOMETRY +++++++++’/ 

A1X’ HOLE’ , 3X’ RADIUS’ , 6X, ’ X’ , 8X, ’ Y’ , 7X, ’ POL. ’ ) 

LPOL=. TRUE. 

LSIZE=. TRUE. 

DIA=2. 0*AI J ( 1 ) 

DO 410 1=1, NHOLE 

WRITE (6, 400) I,AIJ(I),XIJ(I),YIJ(I), PHI JP( I ) 

400 FORMAT! IX, 13, IX, F9. 4, IX, F8. 3, IX, F8. 3, IX, F8. 3) 

AI J ( I )=AIJ ( I )/XLAM 
XI J ( I )=XIJ ( I )/XLAM 
YI J ( I )=YI J ( I )/XLAM 
PHI JP( I )=PHI JP( I )*DTOR 

IF(DABS(AI J ( I )-AI J ( 1 ) ) . GT. 0. 0001 )LSIZE=. FALSE. 

IF (DABS (PH IJP( I )-PHI JP(1 ) ) . GT. 0. 0001 )LPOL=. FALSE. 

410 CONTINUE 

LSAME=. FALSE. 

IF(LSIZE. AND. LPOL)LSAME=. TRUE. 

WRITE(6, 101) 

IF(NUMTE.EQ.O) GO TO 440 

DO 430 1=1 , NUMTE 

WRITE (6, 420) I,MIJ(I),NIJ(I) 

420 FORMAT! IX’ MODE’ 12’ =TE-’I1,I1) 

IF! ( 1+MI J ( I ) ) . GT. 9) GO TO 480 
IF(NIJ(I ). GT. 6) GO TO 480 
IF(MI J ( I ) . LT. 0) GO TO 480 
IF(NI J ( I ) . LT. 1 ) GO TO 480 
430 CONTINUE 
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440 IF(NUMTM.EQ.O) GO TO 470 
DO 460 1=1 , NUMTM 
IP=I+NUMTE 

WRITE (6, 450) IP, MIJP(I),NIJP(I) 

450 FORMAT ( IX’ MODE’ 12’ -TM-’ 11,11) 

IF( ( 1+MIJPf I ) ) . GT. 9) GO TO 480 
IF(NIJP(I).GT. 6) GO TO 480 
IF(MIJP(I).LT.O) GO TO 480 
IF(NIJP( I ) . LT. 1 ) GO TO 480 
460 CONTINUE 
470 WRITE(6, 101 ) 

RETURN 

480 WRITE(6, 490) 

490 FORMAT ( IX’ RANGE OF MODE INDICES EXCEEDED’/ 

A1X’ (M, N ) VALUES ALLOWED: (0,1) THROUGH (8,6)’) 

STOP 

500 WR I TE ( 6 , 5 1 0 ) NHOLE , NMODE , NUMTE , NUMTM 

510 FORMAT ( IX’ NUMBER OF MODES EXCEEDS ARRAY DIMENSIONS’// 
A1X’ NHOLE=’ 13, 5X’ NMODE=’ 13, 5X’ NUMTE=’ 13, 5X’ NUMTM=’ 13) 
STOP 

9999 WRITE(6, 102) 

STOP 

END 
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c. 

c: 

c. 


SUBROUTINE YMAT 


Calculates coefficients CY(II,JJ) of complex admittance matrix [Y] : 


IMPLICIT REAL*8 (A, B, D-H, 0-Z) 

IMPLICIT COMPLEX* 16 (C) 

LOGICAL TEI , TEJ, TMI , TMJ, LPOL, LSIZE, LSAME, LCON , LPLAZ 
DIMENSION Y(2) , YA(2) 

COMMON DIA, MMM, NHOLE, NMODE, NUMTE, NUMTM, F, FCON, ER, LPOL, LSIZE, LSAME 
COMMON /APERS/ IT, MI , NI, MJ, NJ, AI, AJ, TEI , TEJ, TMI , TMJ, PHI , PHIP, R 
COMMON /CONST/ C J , P I , TWOP I , DTOR , RTOD , B IG , SMALL , BMAX , YO , MAX , LAYM 
COMMON /LAYER/ CEP, CUP, NLAY, D(601 ), CE(601 ), CU(601 ), CEUP, LCON, LPLAZ 
COMMON /ARRYS/ MI J (20) , NI J (20) , MI JP(20) , NIJP(20) , 

A CY ( 40 , 40 ) , CA ( 40 , 40 ) , CB ( 40 , 40 ) , IPIV(40), INDX(40,2), 

B AIJ(40),XIJ(40),YIJ(40), PHI JP(40) 

EXTERNAL FCTCWG 
101 FORMAT (IX) 

DO 1000 IH0LE=1 , NHOLE 
DO 1000 JH0LE=1, NHOLE 
DO 1000 IM0DE=1 , NMODE 
DO 1000 JM0DE=1, NMODE 
I I = ( I HOLE- 1 ) * NMODE + I MODE 
JJ= ( JHOLE- 1 ) *NMODE+ JMODE 
IFdl.LE. JJ) GO TO 100 
CY ( 1 1 , J J ) =CY ( J J , II) 

GO TO 1000 
100 CONTINUE 

IFUHOLE.EQ. 1) GO TO 110 
IF(IHOLE.NE. JHOLE) GO TO 110 
IF(. NOT. LSIZE) GO TO 110 
CY (I I , J J ) =CY ( IMODE, JMODE ) 

GO TO 1000 
110 TEI=. TRUE. 

TEJ=. TRUE. 

TMI=. FALSE. 

TMJ=. FALSE. 

IFdMODE.LE. NUMTE) GO TO 120 
TEI=. FALSE. 

TMI=. TRUE. 

120 IF( JMODE. LE. NUMTE) GO TO 130 
TEJ=. FALSE. 

TMJ=. TRUE. 

130 CONTINUE 
IT=0 

IF (TEI. AND. TEJ ) IT=1 
IF(TMI . AND. TMJ ) IT=2 
IF(TEI . AND. TMJ ) IT=3 
IF(TMI. AND. TEJ) IT=4 
AI=AIJ ( IHOLE) 

AJ=AIJ (JHOLE) 
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PHPI=PHIJP( IHOLE) 

PHPJ=PHI JP ( JHOLE ) 

IF(TMI )PHPI=PHPI-0. 5*PI 
I F ( TM J ) PHP J =PHP J -0 . 5*PI 
PHIP=PHPJ-PHPI 
XJI=XIJ( JHOLE)-XIJ ( IHOLE) 
YJI=YIJ(JHOLE)-YIJ( IHOLE) 

R=DSQRT ( XJ I *X J I +YJ I *YJ I ) 

IF ( IHOLE. EQ. JHOLE )R=0. 0 
IF(DABS(R) . LT. 1. 0D-20)R=0. 0 
IF(DABS(PHIP) . LT. 1 . 0D-04)PHIP=0. 0 
PHI=0. 0 

IF(R.LT. 1.0D-04) GO TO 160 
IF (DABS (XJ I ) . LT. 1. 0D-04) GO TO 140 
PH I =DATAN2 ( Y J I , X J I ) -PHP I 
GO TO 150 
140 PHI=0. 5*PI 

IF(YJI. LT. 0. 0)PHI=PHI+PI 
PHI=PHI-PHPI 
150 CONTINUE 
160 IF(TMI) GO TO 170 

IF( . NOT. TEI ) STOP 1111 
MI=MI J ( IMODE) 

NI=NI J ( IMODE) 

GO TO 180 

170 IDEM=IMODE-NUMTE 
MI=MI JP( IDEM) 

NI=NI JP( IDEM) 

180 IF(TMJ) GO TO 190 

IF(. NOT. TEJ) STOP 2222 
MJ=MI J ( JMODE) 

NJ=NI J ( JMODE) 

GO TO 200 

190 JDEM=JMODE-NUMTE 
MJ=MI JP( JDEM) 

NJ=NI JP( JDEM) 

200 CONTINUE 

IF( IHOLE. NE. JHOLE) GO TO 210 
IF(MJ.EQ.MI) GO TO 210 
CY(II, JJ)=(0. 0,0.0) 

GO TO 1000 
210 CONTINUE 

IF( IHOLE. EQ. JHOLE) GO TO 230 
IF ( IHOLE. EQ. 1) GO TO 230 
IF( . NOT. LSIZE) GO TO 230 
IF( . NOT. LPOL) GO TO 230 
IHM1=IH0LE-1 
JHM1= JHOLE- 1 
DO 220 IKH=1 , IHM1 
DO 220 JKH=1, JHM1 
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IF(DABS(AI-AIJ(IKH)).GT. 0.001) GO TO 220 
IF(DABS(AJ-AIJ(JKH)).GT. 0.001) GO TO 220 
TX=DABS(DABS(XJI )-DABS(XI J( IKH)-XIJ( JKH) ) ) 
IF(TX.GT. 0.001) GO TO 220 
TY=DABS ( DABS ( Y J I ) -DABS ( Y I J ( I KH ) - Y I J ( JKH ) ) ) 
IF(TY.GT. 0.001) GO TO 220 

PT=DABS (DABS (PHIP ) -DABS ( PHI JP ( IKH ) -PHI JP ( JKH ) ) ) 

IF(PT.GT. 0.001) GO TO 220 

IK=( IKH-1 )*NMODE+IMODE 

JK- ( JKH-1 ) *NMODE+JMODE 

CY(II, JJ)=CY(IK, JK) 

GO TO 1000 
220 CONTINUE 
230 CONTINUE 
EMI=2. 0 
EMJ=2. 0 

IF(MI. EQ. 0)EMI=1. 0 
IF(MJ. EQ. 0)EMJ=1. 0 

CALL FACS (MI . NI , MJ. NJ, AI , AJ, R, PHI , PHIP, EMI ) 

YTEST=0. 0001 

Y(1)=0. 0 

Y(2)=0. 0 

NDEL=50 

BDEL=0. 25 

AA=0. 0001 

BB=BDEL 

T0L=0. 0001 

240 CALL INTG2 ( AA, BB, FCTCWG, NDEL, TOL, YA, Y ) 

Y ( 1 ) =Y ( 1 ) +YA ( 1 ) 

Y(2)=Y(2)+YA(2) 

IF(NLAY.EQ.O) GO TO 270 
DO 250 N=l, NLAY 
EUREAL=DREAL (CE(N)*CU(N) ) 

IF(EUREAL. LT. 0. 0 ) EUREAL=-EUREAL 
IF ( (BB*BB). LE. EUREAL) GO TO 260 
250 CONTINUE 
GO TO 270 
260 CONTINUE 

IF(DABS(BB-1. 0). LT. 0. 01 )BB=1. 0 
AA=BB 

BB=BB+BDEL 

IF(DABS(AA-1. 0) . LT. 0. 01 )AA=1 . 0001 
I F ( DABS ( BB- 1 . 0 ) . LT . 0 . 0 1 ) BB=0 .9999 
GO TO 240 

270 IF(BB. LT. 2.0) GO TO 260 
BDEL=1. 0 
NBDEL=0 
AA=BB 

BB=BB+BDEL 

280 CALL INTG2 (AA, BB, FCTCWG, NDEL, TOL, YA, Y ) 
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Y( 1 )=Y( 1 )+YA( 1 ) 

Y(2)=Y(2)+YA(2) 

IF(DABS(Y(1 ) ) . LE. YTEST) GO TO 290 
IF(DABS(YA( 1 )/Y ( 1 ) ) . GT. YTEST) GO TO 310 
290 IF(DABS(Y(2)).LE. YTEST) GO TO 300 

IF(DABS(YA(2)/Y(2) ) . GT. YTEST) GO TO 310 
300 NBDEL=NBDEL+ 1 

IF(NBDEL. GE. 5) GO TO 330 
GO TO 320 
310 NBDEL=0 

320 IF(BB.GE.BMAX) GO TO 330 
IF(BB. GE. 20. 0)BDEL=5. 0 
AA=BB 

BB=BB+BDEL 
GO TO 280 

330 CYI J=-YO*DSQRT ( EMI *EMJ ) *DCMPLX ( Y ( 1 ) , Y ( 2 ) ) 

CY(II, JJ)=CYIJ 

WR I TE ( 6 , 340 ) 1 1 , J J , CY ( 1 1 , J J ) , BB , I HOLE , IMODE , JHOLE , JMODE 
340 FORMAT ( IX’ CY(’ 12’ , ’ 12’ )*(’ Ell. 4’ , ' Ell. 4’ ) ’ 3X’ BB=’ F6. 2, 
A4X’ IH=’ 12, 2X’ IM=’ 12, 3X’ JH=’ 12, 2X’ JM=’ 12) 

1000 CONTINUE 

WRITE(6, 101 ) 

RETURN 

END 
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SUBROUTINE SMAT 


c: Calculates coefficients of complex scattering matrix [S] (see eq.l): 


IMPLICIT REAL*8 (A, B, D-H, 0-Z) 

IMPLICIT COMPLEX* 16 (C) 

LOGICAL TEI , TEJ , TMI , TMJ 
LOGICAL LPOL, LSIZE, LSAME 

COMMON DIA, MMM, NHOLE, NMODE, NUMTE, NUMTM, F, FCON, ER, LPOL, LSIZE, LSAME 
COMMON /APERS/ IT, MI , NI, MJ, NJ, AI, AJ, TEI , TEJ, TMI, TMJ, PHI, PHIP, R 
COMMON /ZEROS/ XMNP(9, 6) , XMN(9, 6) 

COMMON /CONST/ CJ, PI , TWOPI , DTOR, RTOD, BIG, SMALL, BMAX, YO, MAX, LAYM 
COMMON /ARRYS/ MI J(20) , NI J (20) , MI JP(20) , NI JP(20) , 

A CY(40, 40) , CA(40, 40) , CB(40, 40) , IPIV(40) , INDX(40, 2) , 

B AIJ(40),XIJ(40),YIJ(40),PHIJP(40) 

DO 130 I HOLE* 1 , NHOLE 
DO 130 JH0LE=1 , NHOLE 
DO 130 I MODE* 1 , NMODE 
DO 130 JM0DE=1, NMODE 
TMI*. FALSE. 

IF( IMODE. GT. NUMTE)TMI=. TRUE. 

I I*( IHOLE-1 ) *NMODE+ IMODE 
JJ* ( JHOLE- 1 ) *NMODE+ JMODE 
AKO I*AI J ( IHOLE ) *TWOPI 
IF(TMI )G0 TO 100 
NI*NIJ( IMODE) 

MIP*MI J ( IM0DE)+1 
FSQI*(XMNP(MIP, NI )/AK0I )**2 
FTSQ=ER-FSQI 
FCON=DABS ( FTSQ ) 

I F ( FCON . LT . SMALL ) FC0N=0 . 0 

IF (FCON. GT. BIG)FCON=BIG 

IF (FTSQ. GE. 0. 0 ) CYMN0=Y0 * DSQRT ( FCON ) 

IF(FTSQ. LT. 0. 0)CYMN0*-CJ*Y0*DSQRT(FCON) 

GO TO 110 

100 I DEM* I MODE-NUMTE 

NI*NIJP( IDEM) 

MIP*MIJP( IDEM)+1 
FSQI*(XMN(MIP, NI )/AKOI )**2 
FTSQ=ER-FSQI 
FCON*DABS ( FTSQ ) 

I F ( FCON . LT . SMALL ) FCON=SMALL 

IF (FCON. GT. BIG)FCON=BIG 

IF (FTSQ. GE. 0. 0)CYMN0*Y0/DSQRT(FC0N) 

IF (FTSQ. LT. 0. 0 )CYMN0=Y0/(-CJ* DSQRT ( FCON ) ) 

110 CONTINUE 

CA(II, JJ)=CY(II, JJ) 

CB(II, JJ)*-CY(II, JJ) 

IFdI.EQ. JJ)CA(II, JJ)=CYMNO+CY(II, JJ) 

IF(II.EQ. JJ)CB(II, JJ)=CYMNO-CY(II, JJ) 
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CY(II, JJ)=CB(II, JJ) 

IF (II. EQ. JJ) WRITE (6, 120) 1 1 , CYMNO, I HOLE, IMODE 
120 FORMAT (lX’YMNC 12’ )=(’ Ell. 4’ , ’ Ell. 4’ )’3X’ IHOLE=’ 12, 3X’ IMODE=’ 12) 
130 CONTINUE 

WRITE(6, 140) 

140 F0RMAT(/10X’ SCATTERING MATRIX’ ) 

CALL CINVCCA, MMM, IPIV, INDX.MAX) 

DO 160 1=1, MMM 
DO 160 J=1 , MMM 
CB(I,J)=(0. 0,0.0) 

DO 150 K=1 , MMM 

CB( I , J)=CB( I , J)+CY ( I, K)*CA(K, J ) 

150 CONTINUE 
160 CONTINUE 

DO 200 1=1, MMM 
DO 200 J=1 , MMM 
XIS0L=-300. 0 
XDX=CDABS ( CB ( I , J ) ) 

IF(XDX. LT. 1. OD-15)GO TO 180 
XISOL=20. *DL0G10(XDX) 

XX1=DIMAG(CB( I, J) ) 

XX2=DREAL ( CB ( I , J ) ) 

PHASE=RT0D*DATAN2 ( XXI , XX2 ) 

WRITEC6, 170)1, J,CB( I, J),XISOL, PHASE 

170 FORMATdX’ S( ’ 12’ , ' 12’ ) = ( ’ Ell. 4’ , ’ Ell . 4’ ) ’ 3XF9. 4’ DB’3XF9.4’ DEG’) 
GO TO 200 

180 WRITEC6, 190)1, J,CB(I,J) 

190 FORMAT! IX’ S( ’ 12’ , ’ 12’ ) = ( ’ Ell . 4* , ’ El 1 . 4’ ) ’ 3X’ BELOW -300 DB’ ) 

200 CONTINUE 
RETURN 
END 
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SUBROUTINE CINV(CA,N, IPIV, INDX, MAX) 

c 

c: Performs matrix inversion of (N by N) complex matrix CA. 

c: On return, CA contains inverted matrix. CA must be dimensioned at 

c: least N by N in calling routine. IPIV and INDX are used internally 
c: by CINV and must be dimensioned at least N and N by 2 in calling 
c: routine. MAX is row dimension of CA and INDX in calling routine. 

c 

IMPLICIT REAL* 8 (A, B, D-H, 0-Z) 

IMPLICIT COMPLEX* 16 (C) 

DIMENSION CA(MAX.N), IPIV(N) , INDX(MAX, 2) 

DATA CO, Cl /(O.ODO.O.ODO), (1. 0D0, 0. 0D0)/ 

DATA BIG, SMALL /I . OD+38, 1 . OD-38/ 

ISCALE=0 
RL=BIG 
RS=SMALL 
CDET=C1 
ADM=1 . 0 
DO 100 J=1 , N 
100 IPIV(J)=0 

DO 220 1=1, N 
AVM=0 . 0 
DO 130 J=1 , N 

IF( IPIV(J) . EQ. 1 ) GO TO 130 
DO 120 K=1 , N 
IF(IPIV(K)-1)110, 120,260 
110 CONTINUE 

AVA=CDABS ( CA ( J , K ) ) 

IF(AVM. GE. AVA) GO TO 120 
IR0W=J 
IC0L=K 
AVM=AVA 
120 CONTINUE 
130 CONTINUE 

IFCAVM.EQ. 0.0) GO TO 250 
IPIV ( IC0L)=IPIV ( IC0L)+1 
IF(IR0W.EQ. ICOL) GO TO 150 
CDET=-CDET 
DO 140 L=1 , N 
CSWAP=CA ( I ROW, L ) 

CA( IROW, L)=CA( ICOL, L) 

CA ( I COL , L ) =CSWAP 
140 CONTINUE 
150 CONTINUE 

INDXU, l)=IROW 
INDX( I , 2)=ICOL 
CPIV=CA( ICOL, ICOL) 

APV=CDABS ( CP I V ) 

IF(APV.EQ.O.O) GO TO 250 
CPIVI=CPIV 


A-18 



ADM=CDABS ( CDET ) 

IF(ADM.LT.RL) GO TO 160 
CDET=CDET/RL 
ADM=CDABS ( CDET ) 

ISCALE=ISCALE+1 
IF(ADH.LT.RL) GO TO 170 
CDET=CDET/RL 
I SCALE= I SCALE+ 1 
GO TO 170 

160 CONTINUE 

IF (ADM. GT. RS) GO TO 170 
CDET=CDET*RL 
ADM=CDABS ( CDET ) 

I SCALE= I SCALE- 1 
IF(ADM.GT.RS) GO TO 170 
CDET=CDET*RL 
I SCALE= I SCALE- 1 

170 CONTINUE 

APV=CDABS ( CP I V I ) 

IF(APV. LT.RL) GO TO 180 
CPIVI=CPIVI/RL 
APV=CDABS ( CP I V I ) 

I SCALE= I SCALE+ 1 
IF(APV. LT.RL) GO TO 190 
CPIVI=CPIVI/RL 
I SCALE= I SCALE+ 1 
GO TO 190 

180 CONTINUE 

IF(APV.GT.RS) GO TO 190 
CPIVI=CPIVI*RL 
APV=CDABS ( CP I V I ) 

I SCALE= I SCALE- 1 
IF(APV.GT.RS) GO TO 190 
CPIVI=CPIVI*RL 
I SCALE= I SCALE- 1 

190 CONTINUE 

CDET=CDET*CPIVI 
CA(ICOL, IC0L)=C1 
DO 200 L=1 , N 

200 CA( ICOL, L)=CA( ICOL, L)/CPIV 
DO 220 Ll=l , N 
IFCL1.EQ. ICOL) GO TO 220 
CSWAP=CA(L1, ICOL) 

CA(L1 , ICOL)=CO 
DO 210 L=l, N 

210 CA (LI , L ) =CA (LI , L ) -CA ( ICOL, L ) *CSWAP 

220 CONTINUE 

DO 240 1=1, N 
L=N+1-I 

IF( INDX(L, 1 ) . EQ. INDX(L, 2) ) GO TO 240 
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IROW=INDX(L, 1) 
ICOL=INDX(L, 2) 

DO 230 K=l, N 
CSWAP=CA(K, IROW) 

CA(K, IROW)=CA(K, ICOL) 
CA(K, ICOL)=CSWAP 
230 CONTINUE 
240 CONTINUE 
GO TO 260 
250 CDET=CO 
ISCALE=0 
260 RETURN 
END 



SUBROUTINE INTG2(AA, BB, FX, NT, TOL, SUM, Y) 


c. 

c: 

c: 

c: 

c: 

c: 

c: 

c: 

c: 

c: 

c: 

c: 

c: 


Integrates 2 real functions (may be real and imaginary parts of 
of a complex function) over limits of AA to BB by using a Simpson 
integration procedure with interval halving until relative error 
TOL is reached on each interval. NT is number of intervals. 

SUM is a 2 element array of the resultant integrals over AA to BB. 

Y is the cumulative values of SUM from previous calls to INTG2. 

Y must be set to zero prior to first call to INTG2 

and is used when integrating FX over wide limits 

by successive calls to INTG2. FX(X,WK) is user : May 1989 

supplied subroutine to evaluate integrands WK. : M.C. Bailey 

FX must be declared EXTERNAL in calling routine. : Hampton, VA 


IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION SUM (2),Y(2),WK(2), SUMA ( 2 ) , SUMB ( 2 ) , ESUM ( 2 ) , BASE ( 2 ) 
DIMENSION A(401 , 2) , B(401 , 2) 

DATA NMAX, NIMAX/400, 2/ 

IF(NT.GT.NMAX) STOP 
DEL=(BB-AA)/NT 
DEL02=0. 5D0*DEL 
DO 100 1=1, NT+1 
X=AA+ ( I - 1 ) *DEL 
CALL FX(X, WK) 

A ( I , 1 ) =WK ( 1 ) 

A ( I , 2 ) =WK ( 2 ) 

100 CONTINUE 

BASE ( 1 ) =DEL02* (A(l , 1 )+A(NT+l, 1 ) ) 

BASE ( 2 ) =DEL02* ( A ( 1 , 2 ) +A ( NT+ 1 , 2 ) ) 

DO 110 1=2, NT 

BASE ( 1 ) =BASE ( 1 ) +DEL* A (1,1) 

BASE ( 2 ) =BASE ( 2 ) +DEL*A (1,2) 

110 CONTINUE 

DEN0M=1 . D-28 
TEST=DABS ( BASE ( 1 ) + Y ( 1 ) ) 

I F ( TEST . GT . DENOM ) DENOM=TEST 
TEST=DABS ( BASE ( 2 ) + Y ( 2 ) ) 

I F ( TEST . GT . DENOM ) DENOM=TEST 

SUM ( 1 ) =0 . 0D0 

SUM(2)=0. 0D0 

DO 300 K=l, NT 

ND=1 


FF=AA+ (K-l ) *DEL 
X=FF +DEL 
CALL FX(X, WK) 

A ( 2 , 1 ) =WK ( 1 ) 

A ( 2 , 2 ) =WK ( 2 ) 

SUMA ( 1 ) =DEL02* ( A ( 1 , 1)+A(2, 1)) 
SUMA ( 2 ) =DEL02* (A(l,2)+A(2,2)) 
200 ND=ND+ND 
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IF(ND.GT.NMAX) GO TO 270 

NB=ND/2 

NA=NB+1 

NP=ND+1 

DELN=DEL/ND 

IF(DELN.LT. l.D-6) GO TO 270 
DO 210 J=l, NB 
X=FF+DELN* ( J+J-l ) 

CALL FX(X.WK) 

B ( J , 1 ) =WK ( 1 ) 

B( J, 2)=WK(2) 

210 CONTINUE 
NAP=NA+1 
NBP=NB+1 
DO 220 J=1 , NA 
JL=NAP-J 
JP=JL+JL-1 
A(JP, 1)=A(JL, 1) 

A( JP, 2)=A( JL, 2) 

220 CONTINUE 

DO 230 J=1 , NB 

JL=NBP-J 

JP=JL+JL 

A( JP, 1 )=B( JL, 1 ) 

A( JP, 2)=B( JL, 2) 

230 CONTINUE 

SUMB( 1 )=0. 5D0* (A( 1, 1 )+A(NP, 1 ) ) 
SUMB(2)=0. 5D0*(A(1, 2)+A(NP, 2) ) 

DO 240 NS=2 , ND 

SUMB( 1 )=SUMB( 1 )+A(NS, 1 ) 

SUMB(2)=SUMB(2)+A(NS, 2) 

240 CONTINUE 

SUMB ( 1 ) =DELN*SUMB ( 1 ) 

SUMB ( 2 ) =DELN*SUMB ( 2 ) 

ESUM ( 1 ) =DABS ( SUMB ( 1 ) -SUMA ( 1 ) ) /DENOM 
ESUM ( 2 ) =DABS ( SUMB ( 2 ) -SUMA ( 2 ) ) /DENOM 
SUMA ( 1 ) =SUMB ( 1 ) 

SUMA(2)=SUMB(2) 

IF ( ESUM ( 1 ) . GT. TOL) GO TO 200 
IF(ESUM(2) . GT. TOL) GO TO 200 
GO TO 280 
270 ND=ND/2 
280 CONTINUE 

SUM ( 1 ) =SUM ( 1 ) +SUMA ( 1 ) 

SUM ( 2 ) =SUM ( 2 ) +SUMA ( 2 ) 

A ( 1 , 1 ) =A ( NP , 1 ) 

A(l, 2)=A(NP, 2) 

300 CONTINUE 
RETURN 
END 
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SUBROUT I NE FCTCWG ( BETA , F I ) 

c 

c: Evaluates integrand (real and Imaginary parts) of equation 2. : 

c: 

IMPLICIT REAL*8 (A, B, D-H, 0-Z) 

IMPLICIT COMPLEX* 16 (C) 

LOGICAL TE I , TE J , TM I , TM J , LCON , LPLAZ 
DIMENSION FI (2) 

COMMON /APERS/ IT, MI , NI , MJ, NJ, AI , AJ, TEI , TEJ, TMI , TMJ, PHI , PHIP, R 
COMMON /LAYER/ CEP , CUP , NLAY ,D(601),CE(601),CU(601), CEUP , LCON , LPLAZ 
CALL MLA YER ( BETA , CW1, CW2) 

CALL UVIJ(BETA, UIJ, VIJ, IT) 

GO TO (10, 20, 30, 40) IT 
GO TO 900 

10 CALL STEI (BETA, SSI, TTI) 

IF (MI. NE. MJ) GO TO 15 

IF(NI.NE.NJ) GO TO 15 

IF(DABS(AI-AJ).GT. l.D-4) GO TO 15 

SSJ=SSI 

TTJ=TTI 

GO TO 50 

15 CALL STEJ (BETA, SSJ, TTJ ) 

GO TO 50 

20 CALL STMI (BETA, SSI, TTI) 

IF(MI.NE.MJ) GO TO 25 

IF ( N I . NE. NJ) GO TO 25 

IF(DABS(AI-AJ).GT. l.D-4) GO TO 25 

SSJ=SSI 

TTJ=TTI 

GO TO 50 

25 CALL STMJ (BETA, SSJ, TTJ) 

GO TO 50 

30 CALL STEI (BETA, SSI, TTI) 

CALL STMJ (BETA, SSJ, TTJ) 

GO TO 50 

40 CALL STMI (BETA, SSI, TTI) 

CALL STEJ ( BETA , SS J , TT J ) 

50 F1=+SSI*SSJ*UI J 

F2=-TTI*TTJ*VI J 
CCW= (F1*CW1+F2*CW2 ) *BETA 
FI ( 1 )=DREAL(CCW) 

FI (2)=DIMAG(CCW) 

RETURN 

900 WRITE(6, 901 ) IT 

901 FORMAT (IX’### ERRORS IN FCTCWG ###’ 

A/1X’ (IT must be between 1 and 4)’ 

B/1X’ IT=’ 15) 

STOP 

END 
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SUBROUTINE MLAYER(BETA, CW1 , CW2 ) 


c: Evaluates complex functions in equations 5 and 6. 


IMPLICIT REAL*8 (A, B, D-H, O-Z) 

IMPLICIT COMPLEX* 16 (C) 

LOGICAL LCON, LPLAZ 

COMMON /LAYER/ CEP, CUP, NLAY, D(601 ), CE(601 ), CU(601 ), CEUP, LCON, LPLAZ 
DATA XKO, BIG, CSMAL/6. 283185307179586D0, 1. OD+28, (1. OD-28, 0. 0D0)/ 
DATA CO, Cl, CJ/(0. 0D0, 0. ODO) , (1. ODO.O. ODO), (0. ODO, 1 . ODO)/ 

CU ( NLAY + 1 ) =CUP 

CE ( NLAY+ 1 ) =CEP 

BSQ=BETA*BETA 

CPF=CO 

CGPOG=CO 

CFP0F=C1 

IF (LCON) GO TO 35 
CPF=C1 

CBSQP=CEUP-BSQ 
I F ( LPLAZ ) CBSQP=C 1 -BSQ 
I F ( DREAL ( CBSQP ) )20, 10, 10 
10 CKZ=XKO*CDSQRT( CBSQP) 

GO TO 30 

20 CKZ=-CJ*XKO*CDSQRT( -CBSQP) 

30 CFPOF =-C J * CKZ 

CGPOG=CFPOF 

35 IF (NLAY. EQ. 0) GO TO 200 

N=NLAY 

40 CARG=CE (N)*CU(N) -BSQ 

IF ( DREAL ( CARG ) . LT. 0. 0) GO TO 50 
CKZ=XKO*CDSQRT ( CARG ) 

GO TO 60 

50 CKZ=-CJ*XKO*CDSQRT ( -CARG ) 

60 CKZD=CKZ*D ( N ) 

A=DREAL ( CKZD ) 

B=DIMAG(CKZD) 

DSINA=DSIN(A) 

DCOSA=DCOS ( A ) 

NSIGN=1 

IF(B. LT. 0. 0)NSIGN=-1 
BP=NSIGN*B 

IF(BP.GT. 1.0D-10) GO TO 70 
CSINKD=DSINA 
CCOSKD=DCOSA 
GO TO 100 

70 IF (BP. LT. 65) GO TO 80 

DC0SHB=0. 5* ( 1 . 6948892D+28 ) 

DS I NHB=NS I GN*DCOSHB 
GO TO 90 

80 EP=DEXP(BP) 
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EM=1 . ODO/EP 
DSINHB=0. 5D0* (EP-EM) 

DCOSHB=0 . 5D0* ( EP+EH ) 
DSINHB=NSIGN*DSINHB 

90 CS INKD-DS INA*DCOSHB+C J * DCOSA*DS I NHB 

CC0SKD=DC0SA*DC0SHB-C J * DS I NA*DS I NHB 
100 CONTINUE 

CCU=CFP0F*(CU(N)/CU(N+1 ) )/CKZ 
CCE=CGP0G*(CE(N)/CE(N+1) )/CKZ 
CFOF 1 ®CPF*CCOSKD-CCU* CS I NKD 
CG0G1 =CCOSKD-CCE*CS I NKD 
CFP0F1=CKZ* (CPF*CSINKD+CCU*CCOSKD ) 
CGP0G1=CKZ* ( CS I NKD+CCE *CCOSKD ) 
CFPOF=CFPOF 1 /CFOF 1 
CGP0G=CGP0G1/CG0G1 
CPF=C1 
N=N-1 

IF(N.GT.O) GO TO 40 
200 CW1»-CJ*XK0*CE(1)/CGP0G 
CW2=CJ*CFP0F/(XK0*CU( 1 ) ) 

RETURN 

END 
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SUBROUTINE FACS(MI I , NI , MJJ, NJ, AI , AJ, R, PHI , PHIP, EMI ) 

IMPLICIT REAL*8 ( A-H, O-Z ) 

COMMON /FACTS/ MI, MJ, MP, MM, MPP, MMP, Mil , MJ1 , MXP, MPI , MPJ, MPPI , MPPJ, 
A EMI , EM2, AIKO, AJKO, RKO, XEI , XEJ, XMI , XMJ, 

B XESI , XESJ , XMSI, XMS J , SF I , SF J , TF I , TF J , TFEI, TFEJ , 

C XDS, XDC, XDSP, XDCP, XDSM, XDCM 

COMMON /ZEROS/ XMNP(9, 6) , XMN(9, 6) 

PI=2. ODO*DASIN( 1 . ODO ) 

TW0PI=2. ODO*P I 

MI=MII 

MJ=MJJ 

MP=MJ+MI 

MM=MJ-MI 

MPP=MP+1 

MMP=MM+1 

MI1=(-1 )**MI 

MJ1=(-1 )**MJ 

MXP=-MM+1 

MPI=MI+1 

MPJ=MJ+1 

MPPI=MPI+1 

MPPJ=MPJ+1 

EM1=EMI-1.0 

EM2=2. O/EMI 

AIKO=AI*TWOPI 

A JKO=A J *TWOP I 

RKO=R*TWOPI 

PHIPMI=MI*PHIP 

PHIPMJ=MJ*PHIP 

XDS=DSIN (PHIPMI ) 

XDC=DCOS(PHIPMI) 

AGMP=MP*PH I -PH I PM J 
AGMM=MM*PH I -PH I PM J 
XDSP=DS I N ( AGMP ) 

XDCP=DCOS ( AGMP ) 

XDSM=DS I N ( AGMM ) 

XDCM=DCOS ( AGMM ) 

XEI=XMNP(MPI , NI ) 

XE J=XMNP (MPJ , N J ) 

XMI=XMN(MPI,NI) 

XM J=XMN ( MPJ , N J ) 

XESI=XEI*XEI 

XESJ =XE J *XE J 

XMSI=XMI*XMI 

XMSJ=XMJ*XMJ 

DEI=XESI-MI*MI 

DEJ=XESJ-MJ*MJ 

FEI=AIKO/DSQRT(DEI ) 

FE J =A JKO/DSQRT ( DEJ ) 

SFI=MI*FEI 
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SFJ=MJ*FEJ 

TFI=XESI*FEI 

TFJ=XESJ*FEJ 

TFEI=-0. 5D0*MI*MI/(XEI*DEI ) 
TFEJ=-0. 5D0*MJ*MJ/(XEJ*DEJ ) 
RETURN 
END 
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SUBROUTINE UVIJ (BETA, UI J, VIJ, IT) 


c: Evaluates functions in table I or table II. 


IMPLICIT REAL*8(A-H, O-Z) 

DIMENSION BJ (21 ) 

COMMON /FACTS/ MI , MJ, MP, MM, MPP, MMP, MI 1 , MJ1 , MXP, MPI , MPJ, MPPI , MPPJ, 
A EMI, EM2, AIKO, AJKO, RKO, XEI, XEJ, XMI, XMJ, 

B XESI, XESJ, XMSI , XMSJ, SFI , SFJ, TFI , TFJ, TFEI , TFEJ, 

C XDS , XDC, XDSP , XDCP , XDSM, XDCM 

UIJ=O.ODO 
VI J=0. ODO 

IF (RKO. GT. 1. 0D-04) GO TO 150 
IF (MI. NE. MJ) RETURN 
GOTO (110, 120, 130, 140)IT 

110 UI J=-EM1*XDC 

VIJ=EM2*XDC 
RETURN 

120 UI J=-EM2*XDC 

RETURN 

130 UIJ=EM1*XDS 
RETURN 

140 UI J=EM1*XDS 
RETURN 

150 ARG=RKO*BETA 

CALL BESJS ( ARG, BJ , MPP ) 

GO TO (160, 190, 220, 250) IT 

160 T1=BJ (MPP ) *XDCP 

IF (MI. GT. MJ) GO TO 170 
T2-MI 1 *BJ (MMP ) *XDCM 
GO TO 180 

170 T2-MJ1*BJ (MXP) *XDCM 

180 UIJ=MJ1* (T1-T2) 

VIJ-MJ1*(T1+T2) 

RETURN 

190 T1=BJ(MPP)*XDCP 

IF (MI. GT.MJ) GO TO 200 
T2=MI 1 *BJ (MMP ) *XDCM 
GO TO 210 

200 T2=M J 1 * B J ( MXP ) * XDCM 

210 UIJ=-MJ1*(T1+T2) 

RETURN 

220 T1=BJ (MPP ) *XDSP 

IF (MI. GT.MJ) GO TO 230 
T2=MI 1 *BJ (MMP ) *XDSM 
GO TO 240 

230 T2=MJ 1 *BJ ( MXP ) *XDSM 

240 UIJ=MJ1* (T1-T2) 

RETURN 

250 T1=BJ(MPP)*XDSP 
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IF(MI.GT.MJ) GO TO 260 
T2-MI 1 *BJ (MMP ) *XDSM 
GO TO 270 

260 T2=MJ 1 *BJ (MXP ) *XDSM 

270 UIJ=MJ1*(T1-T2) 

RETURN 

END 
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SUBROUT I NE STE I ( BETA , SS , TT ) 


c 

c: Evaluates equations 12 and 13 for i-th aperture. 

c: 

IMPLICIT REAL*8 (A-H, 0-Z) 

DIMENSION BJ (21 ) 

COMMON /FACTS/ MI , MJ, MP, MM, MPP, MMP, MI 1 , MJ1 , MXP, MPI , MPJ, MPPI , MPPJ, 
A EMI , EM2, AIKO, AJKO, RKO, XEI , XEJ, XMI , XMJ, 

B XES I , XESJ , XMSI, XMS J , SF I , SF J , TF I , TF J , TFE I , TFEJ , 

C XDS , XDC , XDSP , XDCP , XDSM , XDCM 

ARG=AIKO*BETA 
CALL BES JS ( ARG , B J , MPP I ) 

BJA=BJ (MPI )/ARG 

SS=SFI*BJA 

DN=XESI-ARG*ARG 

IF(DABS(DN).GT. 1.0D-05) GO TO 10 
CALL BESJS(XEI , BJ, MPI ) 

TT=TFEI*BJ (MPI ) 

RETURN 

10 TT=TFI* (MI*BJA-BJ (MPPI ) )/DN 
RETURN 
END 


SUBROUT I NE STM I ( BETA , SS , TT ) 

c 

c: Evaluates equations 14 and 15 for i-th aperture. 

c: 

IMPLICIT REAL*8(A-H, O-Z) 

DIMENSION BJ (21 ) 

COMMON /FACTS/ MI , MJ, MP, MM, MPP, MMP, Mil , MJ1 , MXP, MPI , MPJ, MPPI , MPPJ, 
A EMI , EM2, AIKO, AJKO, RKO, XEI , XEJ, XMI , XMJ, 

B XES I , XESJ, XMSI , XMSJ, SFI, SFJ, TFI, TFJ, TFE I , TFEJ, 

C XDS , XDC , XDSP , XDCP , XDSM , XDCM 

TT=0. ODO 
ARG=AIKO*BETA 
DN=XMSI-ARG*ARG 

IF(DABS(DN).GT. 1.0D-05) GO TO 10 
CALL BESJS(XMI , BJ, MPPI ) 

SS=AIK0*0. 5D0*BJ (MPPI ) 

RETURN 

10 CALL BESJS(ARG, BJ, MPI ) 

SS=AIKO*ARG*BJ (MPI )/DN 

RETURN 

END 
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SUBROUT I NE STE J ( BETA , SS , TT ) 


c: Evaluates equations 12 and 13 for j-th aperture. 


IMPLICIT REAL*8(A-H, 0-Z) 

DIMENSION BJ (21 ) 

COMMON /FACTS/ MI , MJ, MP, MM, MPP, MMP, MI 1 , MJ1 , MXP, MPI , MPJ , MPPI , MPPJ , 
A EMI , EM2, AIKO, AJKO, RKO, XEI , XEJ, XMI , XMJ, 

B XESI, XESJ , XMSI, XMSJ , SF I , SF J , TF I , TF J , TFEI, TFEJ , 

C XDS, XDC, XDSP , XDCP , XDSM , XDCM 

ARG=AJKO*BETA 
CALL BES JS ( ARG , B J , MPPJ ) 

BJA=BJ (MPJ )/ARG 

SS=SFJ*BJA 

DN=XESJ-ARG*ARG 

IF(DABS(DN).GT. 1.0D-05) GO TO 10 
CALL BES JS ( XEJ , B J , MPJ ) 

TT=TFEJ*BJ (MPJ ) 

RETURN 

1 0 TT=TF J * ( M J *B J A-B J ( MPPJ ) ) /DN 

RETURN 
END 


SUBROUT I NE STM J ( BET A , SS , TT ) 


c: Evaluates equations 14 and 15 for j-th aperture. 


IMPLICIT REAL* 8(A-H,0-Z) 

DIMENSION BJ (21 ) 

COMMON /FACTS/ MI , MJ, MP, MM, MPP, MMP, MI 1 , MJ1 , MXP, MPI , MPJ, MPPI , MPPJ, 
A EMI, EM2, AIKO, AJKO, RKO, XEI, XEJ, XMI, XMJ, 

B XESI, XESJ, XMSI, XMSJ, SFI, SFJ, TFI , TFJ, TFEI , TFEJ, 

C XDS , XDC , XDSP , XDCP , XDSM , XDCM 

TT=0. ODO 
ARG=AJKO*BETA 
DN=XMSJ-ARG*ARG 

IF(DABS(DN).GT. 1.0D-05) GO TO 10 
CALL BESJS (XMJ, BJ, MPPJ) 

SS=AJK0*0. 5D0*BJ (MPPJ ) 

RETURN 

10 CALL BESJS (ARG, BJ, MPJ) 

SS=AJKO*ARG*BJ (MPJ )/DN 

RETURN 

END 
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SUBROUT I NE BES JS ( XX , B J , MP ) 

c 

c: Evaluates Bessel function of first kind or orders 0 through MP-1 
c: for argument XX, and stores values in BJ(1) through BJ(MP). 
c: BJ must be dimensioned at least (MP+1) in calling program. 

c: 

IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION B J ( 1 ) , B ( 1 30 ) 

DATA PI.PI02 /3. 14159265358979D0, 1. 57079632679490D0/ 

DATA BIG, SMALL, PI04 /l. OD+38, 1. OD-38, 0. 78539816339745D0/ 

DATA TWOTRD /0. 666666666666667D0/ 

X=DABS(XX) 

BJ ( 1 )=1 . 0D0 

MT=MP-1 

N=MT+1 

IF(X. LT.80. )G0 TO 110 
DO 100 1=1, N 

100 B J ( I ) =DSQRT ( 2 . 0D0/ ( P I *X ) ) *DC0S ( X-P I 04-P I 02* ( I - 1 ) ) 

GO TO 310 
110 CONTINUE 

DO 120 1=2, N 
120 BJ ( I )=0. 0D0 

IF(X. LT. 1. D-06)RETURN 
IF(X-15. ) 130, 130, 140 
130 NTEST=20. 0D0+10. 0D0*X-X**TW0TRD 
GO TO 150 

140 NTEST=90. 0D0+0. 5D0*X 
150 IF(MT-NTEST)170, 160, 160 
160 N=NTEST-1 
GO TO 180 
170 N=MT 
180 BPREV=0. 0 
N1=N+1 
F=2. ODO/X 
D=1 . 0D-06 

IF(X-5. ) 190, 200, 200 
190 MA=X+6 . 0D0 
GO TO 210 

200 MA=1 . 4D0*X+60. ODO/X 

210 IFIXX=X 

MB=N+IFIXX/ 4+2 
MZERO=MAXO (MA, MB ) 

MMAX=NTEST 

DO 280 M=MZERO, MM AX, 3 

FM1=SMALL 

FM=0 . 0D0 

ALPHA=0 . 0D0 

IF (M- (M/2 ) *2 ) 230 , 220, 230 
220 JT=-1 

GO TO 240 
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230 JT=1 
240 M2=M-2 

DO 260 K=1,M2 
MK=M-K 

B(MK)=F*MK*FM1-FM 
IF(B(MK)-BIG)250, 310, 310 
250 FM=FM1 

FM1=B(MK) 

JT=-JT 

S=1+JT 

260 ALPHA=ALPHA+B(MK)*S 
B ( 1 )=F*FM1-FM 
ALPHA=ALPHA+B ( 1 ) 

BTEST=B(N1 ) /ALPHA 

I F ( DABS ( BTEST-BPREV ) -DABS ( D*BTEST ) )290, 290, 270 
270 BPREV=BTEST 
280 CONTINUE 
290 DO 300 1=1, N1 
300 B J ( I ) =B ( I ) /ALPHA 

310 IF (XX. LT. 0. ODO)GO TO 320 
RETURN 
320 N=MT+1 

DO 330 1=1, N 
NN=I-1 

330 BJ ( I )=BJ (I )*(-l. 0)**NN 
RETURN 
END 
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